library(fixest)
library(modelsummary)
library(spmoran)
library(REndo)
library(WeightIt)
library(DirectEffects)
library(vtable)
library(spdep)
library(panelView)
library(ggplot2)
library(ggpubr)
library(viridis)
library(patchwork)
library(cragg)
library(mediation)
library(sensemakr)
library(robomit)
library(spatialreg)
library(rdrobust)

options(scipen = 999) ##scientific notation (no)

data <- readRDS("data_final_merged.rds")  %>% dplyr::select(NUTS_ID, Judet, Year, Month, p_ottoman, share, count, ps_1699, lon, lat, tavg, lograin, logelev, logslope, wheat_rainfed, log_large_rivers_distance,log_sea_lines_distance,route_density_2, vienna_distance, istanbul_distance,geometry)

# Perform regressions with Conley standard errors
model1 <- feols(share ~ p_ottoman , data = data, vcov = 'conley')
model2 <- feols(share ~ p_ottoman + count , data = data, vcov = 'conley')
model3 <- feols(share ~ p_ottoman + count| Year + Month , data = data, vcov = 'conley')
model4 <- feols(share ~ p_ottoman + count + lat + lon + lat*lon | Year + Month , data = data, vcov = 'conley')
model5 <- feols(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin + log_large_rivers_distance + log_sea_lines_distance| Year + Month,  data = data, vcov = 'conley')

msummary(list(model1, model2, model3,model4,model5), stars = c('.' = .1,'*' = .05, '**' = .01, '***' = .001), fmt = "%.5f")

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

p_values <- list(
  model1 = summary(model1)$coeftable[, "Pr(>|t|)"],
  model2 = summary(model2)$coeftable[, "Pr(>|t|)"],
  model3 = summary(model3)$coeftable[, "Pr(>|t|)"],
  model4 = summary(model4)$coeftable[, "Pr(>|t|)"],
  model5 = summary(model5)$coeftable[, "Pr(>|t|)"]
)

lapply(p_values, function(p_vals) {
  -log2(p_vals)
})

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

data <- data %>%
  mutate(unique_id = paste0(NUTS_ID, "_", Year, "_", Month))

unique_data <- data %>%
  dplyr::distinct(unique_id, .keep_all = TRUE)

coords <- cbind(unique_data$lat, unique_data$lon)
nb <- knn2nb(knearneigh(coords, k = 4))  

nb_filtered <- lapply(1:length(nb), function(i) {
  neighbors <- nb[[i]]
  valid_neighbors <- neighbors[unique_data$NUTS_ID[neighbors] != unique_data$NUTS_ID[i]]
  
  if (length(valid_neighbors) == 0) {
    valid_neighbors <- neighbors[1]  
  }
  return(valid_neighbors)
})

attr(nb_filtered, "region.id") <- attr(nb, "region.id")
class(nb_filtered) <- "nb"

listw_filtered <- nb2listw(nb_filtered, style = "W")
listw<-listw_filtered 

# Extract residuals from each model
residuals1 <- residuals(model1)
residuals2 <- residuals(model2)
residuals3 <- residuals(model3)
residuals4 <- residuals(model4)
residuals5 <- residuals(model5)

# Step 8: Calculate Moran's I for each model's residuals
moran1 <- moran.test(residuals1, listw_filtered)
moran2 <- moran.test(residuals2, listw_filtered)
moran3 <- moran.test(residuals3, listw_filtered)
moran4 <- moran.test(residuals4, listw_filtered)
moran5 <- moran.test(residuals5, listw_filtered)

# Print the Moran's I results
list(Moran_I = list(model1 = moran1, model2 = moran2, model3 = moran3, model4 = moran4, model5 = moran5))

# Step 8: Calculate Geary's C for each model's residuals
geary1 <- geary.test(residuals1, listw_filtered)
geary2 <- geary.test(residuals2, listw_filtered)
geary3 <- geary.test(residuals3, listw_filtered)
geary4 <- geary.test(residuals4, listw_filtered)
geary5 <- geary.test(residuals5, listw_filtered)

# Print the Geary's C results
list(Geary_C = list(model1 = geary1, model2 = geary2, model3 = geary3, model4 = geary4, model5 = geary5))

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

model1 <- lm(share ~ p_ottoman , data = data)
model2 <- lm(share ~ p_ottoman + count , data = data)
model3 <- lm(share ~ p_ottoman + count + Year + Month , data = data)
model4 <- lm(share ~ p_ottoman + count + lat + lon + Year + Month , data = data)
model5 <- lm(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin + log_large_rivers_distance + log_sea_lines_distance + Year + Month,  data = data)

# runs sensemakr for sensitivity analysis
sensitivity2 <- sensemakr(model = model2, treatment = "p_ottoman",benchmark_covariates = "count",kd = 1:3)
sensitivity3 <- sensemakr(model = model3, treatment = "p_ottoman",benchmark_covariates = "count",kd = 1:3)
sensitivity4 <- sensemakr(model = model4, treatment = "p_ottoman",benchmark_covariates = "count",kd = 1:3)
sensitivity5 <- sensemakr(model = model5, treatment = "p_ottoman",benchmark_covariates = "count",kd = 1:3)

summary(sensitivity2)
summary(sensitivity3)
summary(sensitivity4)
summary(sensitivity5)

########################################################################################################################################################################################
########################################################################################################################################################################################
########################################################################################################################################################################################

# Tables A1-A3

# Perform regressions with clustered standard errors (Year-Month)
model1 <- feols(share ~ p_ottoman , data = data, cluster = ~ Year + Month )
model2 <- feols(share ~ p_ottoman + count , data = data, cluster = ~ Year + Month )
model3 <- feols(share ~ p_ottoman + count| Year + Month , data = data, cluster = ~ Year + Month )
model4 <- feols(share ~ p_ottoman + count + lat + lon + lat*lon | Year + Month , data = data, cluster = ~ Year + Month )
model5 <- feols(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin + log_large_rivers_distance + log_sea_lines_distance| Year + Month,  data = data, cluster = ~ Year + Month )

msummary(list(model1, model2, model3,model4,model5), stars = c('.' = .1,'*' = .05, '**' = .01, '***' = .001), fmt = "%.5f", output = "latex")

# Perform regressions with clustered standard errors (NUTS)
model1 <- feols(share ~ p_ottoman , data = data,cluster = ~ NUTS_ID )
model2 <- feols(share ~ p_ottoman + count , data = data, cluster = ~ NUTS_ID )
model3 <- feols(share ~ p_ottoman + count| Year + Month , data = data, cluster = ~ NUTS_ID )
model4 <- feols(share ~ p_ottoman + count + lat + lon + lat*lon | Year + Month , data = data, cluster = ~ NUTS_ID )
model5 <- feols(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin + log_large_rivers_distance + log_sea_lines_distance| Year + Month,  data = data, cluster = ~ NUTS_ID )

msummary(list(model1, model2, model3,model4,model5), stars = c('.' = .1,'*' = .05, '**' = .01, '***' = .001), fmt = "%.5f", output = "latex")

# Perform regressions with DK standard errors
model1 <- feols(share ~ p_ottoman , data = data, vcov = 'DK', panel.id = ~Judet + Year)
model2 <- feols(share ~ p_ottoman + count , data = data,  vcov = 'DK', panel.id = ~Judet + Year)
model3 <- feols(share ~ p_ottoman + count| Year + Month , data = data,  vcov = 'DK', panel.id = ~Judet + Year)
model4 <- feols(share ~ p_ottoman + count + lat + lon + lat*lon | Year + Month , data = data,  vcov = 'DK', panel.id = ~Judet + Year)
model5 <- feols(share ~ p_ottoman + count + ps_1699 + tavg + logelev + logslope + lograin + log_large_rivers_distance + log_sea_lines_distance| Year + Month,  data = data,  vcov = 'DK', panel.id = ~Judet + Year)

msummary(list(model1, model2, model3,model4,model5), stars = c('.' = .1,'*' = .05, '**' = .01, '***' = .001), fmt = "%.5f", output = "latex")

